# Generazione sinusoidale con diverse armoniche
signal <- function(t, pars, rad = FALSE) {
stopifnot(is.data.frame(pars))
with(pars, {
if (!rad) {
phi <- phi/180*pi
f <- 2*pi*f
}
map_dbl(t, \(t) sum( map_vec(seq_along(w) , \(i) w[i]*sin(t*f[i] + phi[i] ))))
})
}Esercitazione 4
Compensazione dinamica
Per la consegna dell’esercitazione si faccia riferimento a questo link. Mentre le librerie sono rimandate nell’icona Librerie.
1 Funzioni di utilità
2 Simulazione uscita sistema di misura
# Parametri del misurando (input del sistema di misura)
pars <- tibble(
w = c(1, 0.3, 0.3),
f = c(1/30, 1/10, 1/4),
phi = c(0, 0, 0)
)
# Generiamo il misurando u
Ts <- 0.01 # periodo di campionamento
fs <- 1/Ts # freq. di campionamento
u <- tibble(
t = 0:10000 * Ts,
u = signal(t, pars),
un = u + rnorm(length(t), 0, pars$w[1]/10)
)# Definiamo la Funzione di Trasferimento dello strumento
num <- c(5, 1) # 5s + 1
den <- c(1, 2, 1) # s^2 + 2s + 1
H <- tf(num, den)
print(H)
y1:
5 s^1 + 1
- - - - - - - - - -
s^2 + 2 s + 1
Transfer Function: Continuous time model

# Simulazione dell'uscita nel tempo
tmp <- u %>% mutate(
out = lsim(H, un, t)$y[1,],
ideal = lsim(H, u, t)$y[1,]
) %>% mutate(n = row_number(), .before = t)
tmp %>% plot_ly(x = ~t) %>%
add_lines(y = ~un, name = "Input") %>%
add_lines(y = ~ideal, name = "Output ideale") %>%
add_lines(y = ~out, name = "Output") %>%
layout(xaxis=list(title="Tempo (s)"), yaxis=list(title=""), title = "Simulazione dell'uscita")Notare che lsim calcola la risposta nel tempo di un sistema lineare.
2.1 Spettro dei segnali misurati
Ta <- max(u$t)
un.fft <- compute_fft(tmp, Ta, un)
ideal.fft <- compute_fft(tmp, Ta, ideal)
out.fft <- compute_fft(tmp, Ta, out)
bind_rows(
un.fft %>% mutate(caso = "Input"),
ideal.fft %>% mutate(caso = "Output ideale"),
out.fft %>% mutate(caso = "Output")
) %>% dplyr::filter(f < 1) %>%
ggplot(aes(x = f, y = mod)) +
geom_spoke(aes(angle = pi/2, y=0, radius = mod), linewidth = 0.2, colour = "#F8766D") +
geom_point(colour = "#F8766D") +
facet_wrap(~ caso, ncol = 1, scales = "free_x") +
labs(
x = "Frequenza (Hz)",
y = "Intensità"
)
NB: non si vede, ma il rumore bianco è presente. Provare ad ingrandire la dev. std. del rumore addizionato a un.
3 Stima del misurando a partire dal segnale di uscita
# Segnale in uscita con rumore
yI <- lsim(H, u$u, u$t)$y[1,]
yn <- yI + rnorm(length(u$t), 0, pars$w[1]/10)
# Ricaviamo l'inversa della Funzione di Trasferimento:
Hinv <- tf(den, num)[1] "TFCHK: Transfer function may not be proper and may lead to errors. Num > Den"
print(Hinv)
y1:
s^2 + 2 s + 1
- - - - - - - - - -
5 s^1 + 1
Transfer Function: Continuous time model

# Simulazione dell'ingresso
#input <- lsim(Hinv, yn, y$t, x0 = rep(0, length(pole(Hinv))))Error in tf2ss(num, den) :
TF2SS: The order of the Numerator should be equal or lesser than the Denominator.
Implementazione della compensazione dinamica.
Il problema risiede nel fatto che la funzione lsim è progettata per calcolare la risposta nel tempo di un sistema LTI proprio (\(deg(\text{NUM}) \leq deg(\text{DEN})\)). Nel caso dell’inversa \(H^{-1}(j\omega)\), è come se si stimasse la risposta di un sistema improprio, per cui lsim fallisce.
Risolvo elaborando tutto nel dom della freq e poi antitrasformo.
# filtro l'uscita per ridurre il rumore
y <- u %>% select(t) %>% mutate(n = row_number(t), .before = t) %>%
mutate(
y_id = yI,
y_n = yn
)
compute_fft(y, Ta, yn) %>% dplyr::filter(f < 2.5) %>%
ggplot(aes(x = f, y = mod)) +
geom_spoke(aes(angle = pi/2, y=0, radius = mod), linewidth = 0.2, colour = "#56B4E9") +
geom_point(colour = "#56B4E9") +
labs(x = "Frequenza (Hz)", y = "Intensità", title = "Spettro dell'uscita con rumore")
# filtro a 1Hz
fc <- 1
flt <- butter(2, fc/(fs/2))
ggbodeplot_digital(flt, fs, fmin = 1e-1, fmax = 1e1) + geom_vline(xintercept = fc, color="red", linetype = 2)
y <- y %>% mutate(y_flt = signal::filtfilt(flt, yn))
# per padding automatico di filtfilt ho una distorsione, perciò taglio il segnale di uscita
y <- y %>% dplyr::filter(t <= 99)
# Nuovo tempo di campionamento
Ta <- max(y$t)# Windowing perché segnale complessivo non periodico -> FFT non funziona
N <- length(y$t)
y <- y %>% mutate(
win = signal::hanning(N),
y_win = y_flt * win
)
y %>% pivot_longer(cols = c("y_win", "y_flt")) %>%
ggplot(aes(x = t, y = value, colour = name)) +
geom_line() +
geom_line(aes(y=win), color = "red", linetype = 22) +
geom_line(aes(y=-win), color = "red", linetype = 22) +
scale_color_discrete(labels = c(y_win = "Windowing", y_flt = "Filtrata")) +
labs(x = "Tempo (s)", y="Output", colour = "Segnale")
[ChatGPT] Per chiarezza e per evitare errori numerici, viene calcolata esplicitamente la funzione di risposta in frequenza, anziche usare la freqresp.
# Trasformata di Fourier (FFT)
y <- y %>% mutate(
Y =fft(y_flt),
Y_win = fft(y_win)
)
# Asse frequenze
freq <- tibble(
k = 0:(N-1),
f = ifelse(k <= N/2, k, k-N) * fs/N, # FFT simmetrica nell'asse delle freq.
w = 2*pi*f
)
# Calcolo della funz. di risposta in frequenza
Hiw <- freqresp(H, freq$w)# Inverto la caratteristica dinamica e anti-trasformo
u_comp <- tibble(
t = y$t,
U = y$Y / Hiw,
U_win = y$Y_win / Hiw,
u = Re(fft(U, inverse = TRUE)) / N,
u_win = Re(fft(U_win, inverse = TRUE)) / N # Re() serve a trascurare dei residui di calcolo nella parte Im
# la fft(inverse = TRUE) ritorna l'inversa NON-normalizzata, perciò si divide per la lunghezza del campione
)u_comp %>% plot_ly(x = ~t) %>% #add_lines(y = y$y_n, name = "Uscita") %>%
add_lines(y = u$un[u$t <= Ta], name = "Reale") %>%
add_lines(y = u$u[u$t <= Ta], name = "Ideale") %>%
add_lines(y = ~u, name = "Compensato") %>%
add_lines(y = ~u_win, name = "Windowing")3.1 Provo con padding simmetri a zero [ChatGPT]
# Ampiezza padding
M <- 4*N
pad_left <- floor((M - N)/2)
pad_right <- ceiling((M - N)/2)
t_start <- y$t[1] - pad_left * Ts
y_ext <- tibble(
t = seq(from = t_start, by = Ts, length.out = M),
y_pad = c(
rep(0, pad_left),
y$y_flt,
rep(0, pad_right)
)
)
# Verifica che il segnale sia centrato: deve restituire TRUE
all.equal(
y_ext$y_pad[(pad_left+1):(pad_left+N)],
y$y_flt
)[1] TRUE
y_ext %>%
ggplot(aes(t, y_pad)) +
geom_line() +
geom_line(data = y, aes(x = t, y = y_n), color = "red")
freq <- tibble(
k = 0:(M-1),
f = ifelse(k <= M/2, k, k-M) * fs/M, # FFT simmetrica nell'asse delle freq.
w = 2*pi*f
)
Hiw <- freqresp(H, freq$w)
y_ext <- y_ext %>% mutate(
Y_pad = fft(y_pad),
U_pad = Y_pad/Hiw,
u_pad = Re(fft(U_pad, inverse = TRUE)) / M
)
u_comp <- u_comp %>% mutate(
u_pad = y_ext$u_pad[(pad_left+1):(pad_left+N)]
)u_comp %>% plot_ly(x = ~t) %>% #add_lines(y = y$y_n, name = "Uscita") %>%
add_lines(y = u$un[u$t <= Ta], name = "Reale") %>%
add_lines(y = u$u[u$t <= Ta], name = "Ideale") %>%
add_lines(y = ~u, name = "Compensato") %>%
add_lines(y = ~u_win, name = "Windowing") %>%
add_lines(y = ~u_pad, name = "Padded")Applico un filtro sul segnale padded:
compute_fft(u_comp, Ta, u_pad) %>% dplyr::filter(f < 2) %>%
ggplot(aes(x = f, y = mod)) +
geom_spoke(aes(angle = pi/2, y=0, radius = mod), linewidth = 0.2, colour = "#56B4E9") +
geom_point(colour = "#56B4E9") +
labs(x = "Frequenza (Hz)", y = "Intensità", title = "Spettro dell'ingresso compensato")
fc <- 2
flt <- butter(2, 2*fc/fs)
u_comp <- u_comp %>% mutate(
u_flt = signal::filtfilt(flt, u_pad)
)
u_comp %>% plot_ly(x = ~t) %>%
add_lines(y = u$un[u$t <= Ta], name = "Reale") %>%
add_lines(y = u$u[u$t <= Ta], name = "Ideale") %>%
add_lines(y = ~u_pad, name = "Padded") %>%
add_lines(y = ~u_flt, name = "Filtered")ora accetto il risultato